% y代表原序列
[y,Fs]=audioread('test.wav');
y=y(:,1);
T=1/Fs;
% sound(y,Fs)
t=(1:length(y))*T;
subplot(3,2,1)
plot(t,y)
ylim([-0.6 0.6])
xlabel("时间/s");
ylabel("幅度");
title("原始语音信号时域图");
Y=fft(y);
subplot(3,2,2)
plot((1:length(Y))*2/length(Y),abs(Y))
xlabel("频率");
ylabel("幅度");
title("原始语音信号频域图");
ylim([0 600])
% 加入高斯白噪声，即加噪后的序列y1
y1=y+0.2*rand(size(y))-0.1;
%sound(y1,Fs)
% audiowrite("加噪.wav",y1,Fs);
subplot(3,2,3)
plot(t,y1)
ylim([-0.6 0.6])
xlabel("时间/s");
ylabel("幅度");
title("加噪音语音信号时域图");
Y1=fft(y1);
subplot(3,2,4)
plot((1:length(Y1))*2/length(Y1),abs(Y1))
ylim([0 600])
xlabel("频率");
ylabel("幅度");
title("加噪音语音信号频域图");
% 降噪处理
% 巴特沃斯低通滤波器
%通带截止频率
Wp=0.20;
% 阻带截止频率
Ws=0.55;
Rp=3;Rs=30;
% Wp=2*pi*fp;Ws=2*pi*fs;

[N,Wc]=buttord(Wp,Ws,Rp,Rs);
[b,a]=butter(N,Wc,"low");


y2 = filter(b,a,y1);
% sound(y2, Fs);

% 去噪音语音信号时域图
subplot(3,2,5);
plot(t,y2);
ylim([-0.6 0.6]);
xlabel("时间/s");
ylabel("幅度");
title("巴特沃斯去噪音语音信号时域图");

Y2 = fft(y2);
% 去噪音语音信号频域图
subplot(3,2,6);
plot((1:length(Y2))*2/length(Y2),abs(Y2));
ylim([0 600]);
xlabel("频率");
ylabel("幅度");
title("巴特沃斯去噪音语音信号频域图");

% 计算滤波前的信噪比
% 原始信号功率

psignal=P(y);
% 噪声功率
pnoise=P(y1-y);
% 信噪比
SNR_Original=10*log10(psignal/pnoise);
fprintf('滤波前信噪比为：%.3f\n',SNR_Original);
% disp(SNR_Original)

% 计算滤波后的信噪比
% 降噪后的功率
psignal_after=P(y);
% 噪声功率
pnoise=P(y2-y);
SNR_After=10*log10(psignal_after/pnoise);
% disp(SNR_After)
fprintf('滤波后信噪比为：%.3f\n',SNR_After);

% 新的图形窗口绘制滤波器频率响应图
figure
plot(w/pi,20*log10(abs(H)));
xlabel('频率');
ylabel('幅度 (dB)');
title('巴特沃斯低通滤波器频率响应');


